///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using ILNumerics;
using ILNumerics.BuiltInFunctions;

// Linear Equations functions
// ==========================
//
// chol - Cholesky factorization 
// cholinc - Sparse incomplete Cholesky and Cholesky-Infinity factorizations 
// cond - Condition number with respect to inversion 
// condest - 1-norm condition number estimate 
// funm - Evaluate general matrix function 
// ilu - Sparse incomplete LU factorization 
// inv - Matrix inverse 
// ldl - Block LDL' factorization for Hermitian indefinite matrices 
// linsolve - Solve linear system of equations 
// lscov - Least-squares solution in presence of known covariance 
// lsqnonneg - Solve nonnegative least-squares constraints problem 
// lu - LU matrix factorization 
// luinc - Sparse incomplete LU factorization 
// pinv - Moore-Penrose pseudoinverse of matrix 
// qr - Orthogonal-triangular decomposition 
// rcond - Matrix reciprocal condition number estimate 

namespace MathLib
{
    /// <summary>
    /// Linear Equations functions.
    /// </summary>
    public partial class MLMath
    {
        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<complex> inv(ILArray<complex> a1)
        {
            ILArray<complex> outArray = new ILArray<complex>(a1.Dimensions.ToIntArray());

            MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Complex.DenseMatrix(a1.Dimensions[1], a1.Dimensions[0]);
            for (int r = 0; r < a1.Dimensions[1]; r++)
                for (int c = 0; c < a1.Dimensions[0]; c++)
                    a[r, c] = new System.Numerics.Complex(
                        a1.InternalArray4Experts[r * a1.Dimensions[0] + c].real,
                        a1.InternalArray4Experts[r * a1.Dimensions[0] + c].imag);

            MathNet.Numerics.LinearAlgebra.Generic.Matrix<System.Numerics.Complex> minv = a.Inverse();

            for (int r = 0; r < outArray.Dimensions[1]; r++)
                for (int c = 0; c < outArray.Dimensions[0]; c++)
                    outArray.InternalArray4Experts[r * outArray.Dimensions[0] + c] = new complex(
                        minv[r, c].Real,
                        minv[r, c].Imaginary);

            return outArray;
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<fcomplex> inv(ILArray<fcomplex> a1)
        {
            ILArray<fcomplex> outArray = new ILArray<fcomplex>(a1.Dimensions.ToIntArray());

            MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Complex32.DenseMatrix(a1.Dimensions[1], a1.Dimensions[0]);
            for (int r = 0; r < a1.Dimensions[1]; r++)
                for (int c = 0; c < a1.Dimensions[0]; c++)
                    a[r, c] = new MathNet.Numerics.Complex32(
                        a1.InternalArray4Experts[r * a1.Dimensions[0] + c].real,
                        a1.InternalArray4Experts[r * a1.Dimensions[0] + c].imag);

            MathNet.Numerics.LinearAlgebra.Generic.Matrix<MathNet.Numerics.Complex32> minv = a.Inverse();

            for (int r = 0; r < outArray.Dimensions[1]; r++)
                for (int c = 0; c < outArray.Dimensions[0]; c++)
                    outArray.InternalArray4Experts[r * outArray.Dimensions[0] + c] = new fcomplex(
                        minv[r, c].Real,
                        minv[r, c].Imaginary);

            return outArray;
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<double> a1)
        {
            ILArray<double> outArray = new ILArray<double>(a1.Dimensions.ToIntArray());

            MathNet.Numerics.LinearAlgebra.Double.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(a1.Dimensions[1], a1.Dimensions[0]);
            for (int r = 0; r < a1.Dimensions[1]; r++)
                for (int c = 0; c < a1.Dimensions[0]; c++)
                    a[r, c] = a1.InternalArray4Experts[r * a1.Dimensions[0] + c];

            MathNet.Numerics.LinearAlgebra.Generic.Matrix<double> minv = a.Inverse();

            for (int r = 0; r < outArray.Dimensions[1]; r++)
                for (int c = 0; c < outArray.Dimensions[0]; c++)
                    outArray.InternalArray4Experts[r * outArray.Dimensions[0] + c] = minv[r, c];

            return outArray;
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<float> inv(ILArray<float> a1)
        {
            ILArray<float> outArray = new ILArray<float>(a1.Dimensions.ToIntArray());

            MathNet.Numerics.LinearAlgebra.Single.DenseMatrix a = new MathNet.Numerics.LinearAlgebra.Single.DenseMatrix(a1.Dimensions[1], a1.Dimensions[0]);
            for (int r = 0; r < a1.Dimensions[1]; r++)
                for (int c = 0; c < a1.Dimensions[0]; c++)
                    a[r, c] = a1.InternalArray4Experts[r * a1.Dimensions[0] + c];

            MathNet.Numerics.LinearAlgebra.Generic.Matrix<float> minv = a.Inverse();

            for (int r = 0; r < outArray.Dimensions[1]; r++)
                for (int c = 0; c < outArray.Dimensions[0]; c++)
                    outArray.InternalArray4Experts[r * outArray.Dimensions[0] + c] = minv[r, c];

            return outArray;
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<long> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<ulong> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<int> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<uint> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<short> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<ushort> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<char> a1)
        {
            return inv(ILMath.todouble(a1));
        }

        /// <summary>
        /// Computes the inverse of a square matrix
        /// </summary>
        /// <param name="a1">Square matrix to compute inverse for</param>
        /// <returns>Inverse of square matrix</returns>
        public static ILArray<double> inv(ILArray<byte> a1)
        {
            return inv(ILMath.todouble(a1));
        }
    }
}